home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / powers.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  10.1 KB  |  305 lines
  1. ;; Maxima code for extracting powers, finding leading and trailing 
  2. ;; coefficients, and finding the degree of polynomials.
  3.  
  4. ;; Author Barton Willis, University of Nebraska at Kearney (aka UNK)
  5. ;; December 2001, December 2002
  6.  
  7. ;; License: GPL
  8. ;; The user of this code assumes all risk for its use. It has no warranty.
  9. ;; If you don't know the meaning of "no warranty," don't use this code. :)
  10.  
  11. (in-package "MAXIMA")
  12. ($put '$powers 1 '$version)
  13.  
  14. ;; Acknowledgement: Dan Stanger helped find and correct bugs.  He
  15. ;; also wrote user documentation and a test routine. 
  16.  
  17. ;; emptyp(e) returns true iff e is an empty Maxima list. This function 
  18. ;; is identical to is(e = []).  This function isn't used by any functions 
  19. ;; in powers; it could be expunged.
  20.  
  21. (defun $emptyp (e)
  22.   (like e '((mlist))))
  23.  
  24. ;; posintp(x) returns true iff x is a positive integer or if x has been declared 
  25. ;; to be an integer and has been assumed to be greater than zero. Thus
  26.  
  27. ;;  (C1) declare(m, integer)$
  28. ;;  (C2) assume(m > 0)$
  29. ;;  (C3) posintp(m);
  30. ;;  (D3)        TRUE 
  31. ;; posintp isn't used by any functions in powers; it could be expunged.
  32.  
  33. (defun $posintp (x)
  34.   (and (or ($integerp x) ($featurep x '$integer)) (mgrp x 0)))
  35.  
  36. ;; Set ratfac to nil, return rat(e,x), and reset ratfac to
  37. ;; its previous value.
  38.  
  39. (defun myrat (e &rest x)
  40.   (let ((save-ratfac $ratfac))
  41.     (setq $ratfac nil)
  42.     (unwind-protect
  43.     (apply '$rat `(,e ,@x))
  44.       (setq $ratfac save-ratfac))))
  45.             
  46. ;; If x list a Maxima list of symbols, return true iff the expression p
  47. ;; doesn't depend on any member of x.  If x is a symbol, return true
  48. ;; iff p doesn't depend on x.  This function is similar to $freeof, but
  49. ;; it maybe somewhat more efficient when p is a CRE expression. Finally,
  50. ;; if x (any member of x when x is a Maxima list) isn't a symbol signal 
  51. ;; an error.
  52.  
  53. (defun $ratfreeof (x p)
  54.   (setq x (require-list-of-symbols x "$ratfreeof" 2))
  55.   (let ((p-vars (cdr ($showratvars p))))
  56.     (cond ((every #'(lambda (z) (or ($symbolp z) ($subvarp z))) p-vars)
  57.        (every #'(lambda (z) (null (member z p-vars :test #'like))) x))
  58.       (t 
  59.        (setq p ($totaldisrep p))
  60.        (every #'(lambda(z) ($freeof ($totaldisrep z) p)) x)))))
  61.  
  62. ;; variablep(e) evaluates to true if and only if e is a non-constant symbol 
  63. ;; or a subscripted symbol. Because symbolp(pi) evaluates to true, we need to 
  64. ;; check whether cd mae is constant.
  65.  
  66. (defun $variablep (e)
  67.   (and (or ($symbolp e) ($subvarp e)) (not ($constantp e))))
  68.  
  69. ;; ordinal_string(i) returns the ordinal name of the integer i.  When 
  70. ;; i > 10, i < 1, or i isn't an integer, give up and return i-th.
  71.  
  72. (defun $ordinal_string (i)
  73.   (case i
  74.     (1 "first")
  75.     (2 "second")
  76.     (3 "third")
  77.     (4 "fourth")
  78.     (5 "fifth")
  79.     (6 "sixth")
  80.     (7 "seventh")
  81.     (8 "eighth")
  82.     (9 "ninth")
  83.     (10 "tenth")
  84.     (otherwise (format nil "~A-th" 
  85.                (string-left-trim "&" (mfuncall '$string i))))))
  86.     
  87. ;; If variablep(v) evaluates to false, signal an error saying that
  88. ;; the i-th argument of the function f requires a symbol; otherwise, 
  89. ;; return true.
  90.  
  91. (defun require-symbol (v f i)
  92.   (if (not ($variablep v))
  93.       (merror "The ~A argument of ~:M must be a symbol, instead found ~:M" 
  94.           ($ordinal_string i) f v) t))
  95.  
  96. ;; If v is a Maxima list and each element of v is a symbol, return the
  97. ;; cdr of v.  When v isn't a list, but is a symbol, return the Lisp list
  98. ;; (list v). Otherwise signal an error saying that the i-th argument of the
  99. ;; function f requires a symbol or a list of symbols.
  100.  
  101. (defun require-list-of-symbols (v f i)
  102.   (let ((x))
  103.     (if ($listp v) (setq x (cdr v)) (setq x (list v)))
  104.     (if (every #'$variablep x) x
  105.       (merror "The ~A argument of ~:M must be a symbol or a list of symbols, instead found ~:M" ($ordinal_string i) f v))))
  106.       
  107. (defun require-poly (p v f i)
  108.   (setq p (myrat p v))
  109.   (if ($polynomialp p v) p
  110.     (merror "The ~A argument of ~:M requires a polynomial, instead found ~:M" ($ordinal_string i) f p)))
  111.  
  112. (defun require-nonlist (e f i)
  113.   (if ($listp e)
  114.       (merror "The ~A argument of ~:M requires a nonlist, instead found ~:M" 
  115.           ($ordinal_string i) f e)))
  116.  
  117. ;; Return a Maxima list of the non-constant rat variables in e.
  118.    
  119. (defun non-constant-ratvars (e)
  120.   (let ((v (cdr ($showratvars e)))
  121.     (acc))
  122.     (dolist (vi v `((mlist simp) ,@acc))
  123.       (if (not ($constantp vi)) (push vi acc)))))
  124.  
  125. ;; If e is a list, map $powers over the list.  If e is a sum of powers 
  126. ;; of powers of x, return a list of the exponents.
  127.  
  128. (defun $powers (e x)
  129.   (require-symbol x "$powers" 2)
  130.   (cond (($listp e)
  131.      (cons '(mlist simp) (mapcar #'(lambda (p) ($powers p x)) (cdr e))))
  132.     (t
  133.      (setq e (require-poly (myrat e x) x "$powers" 1))
  134.      (cond (($ratfreeof x e)
  135.         `((mlist simp) 0))
  136.            (t
  137.         (cons '(mlist simp) (odds (cadr e) 0)))))))
  138.            
  139. ;; odds is defined in mactex.lisp. Here is its definition.
  140.  
  141. (defun odds (n c)
  142.   ;; if c = 1, get the odd terms  (first, third...)
  143.   (cond ((null n) nil)
  144.     ((= c 1) (cons (car n) (odds (cdr n) 0)))
  145.     ((= c 0) (odds (cdr n) 1))))
  146.  
  147. ;; Return the highest power of the polynomial e in the variable x.
  148.  
  149. (defun $hipower (e x)
  150.   (require-symbol x "$hipower" 2)
  151.   (setq e (require-poly (myrat e x) x "$hipower" 1))
  152.   (if (or ($constantp e) ($ratfreeof x e)) 0 (cadadr e)))
  153.  
  154. ;; Return the lowest power of the polynomial e in the variable x.
  155.  
  156. (defun $lowpower (e x)
  157.   (require-symbol x "$lowpower" 2)
  158.   (setq e (require-poly (myrat e x) x "$lowpower" 1))
  159.   (if (or ($constantp e) ($ratfreeof x e)) 0 (nth 1 (reverse (cadr e)))))
  160.  
  161. ;; Flatten a Maxima list.
  162.  
  163. (defun flatten-list (e)
  164.   (cond (($listp e)
  165.      (let ((acc))
  166.        (dolist (ei (cdr e) (cons '(mlist simp) (nreverse acc)))
  167.          (setq acc (if ($listp ei) (nconc (cdr (flatten-list ei)) acc)
  168.              (cons ei acc))))))
  169.     (t e))) 
  170.  
  171. ;; If e is a sum of powers of x, return a list of the coefficients
  172. ;; of powers of x.  When e isn't a sum of powers, return false.  This 
  173. ;; function is based on a Macsyma function written by A.D. Kennedy and 
  174. ;; referenced in "Mathematics and System Reference Manual," 16th edition, 
  175. ;; 1996.
  176.  
  177. (defun $allcoeffs (e x)
  178.   (flatten-list (allcoeffs e x)))
  179.  
  180. (defun allcoeffs (e x)
  181.   (cond (($listp e)
  182.      (cons '(mlist simp) (mapcar #'(lambda (s) (allcoeffs s x)) (cdr e))))
  183.     (($listp x)
  184.      (cond ((= 0 ($length x)) e)
  185.            (t (allcoeffs (allcoeffs e ($first x)) ($rest x)))))
  186.     (t 
  187.      (require-symbol x "$allcoeffs" 2)
  188.      (setq e (myrat e x))
  189.      (let ((p ($powers e x)))
  190.        (cons '(mlist simp) 
  191.          (mapcar #'(lambda (n) ($ratcoef e x n)) (cdr p)))))))
  192.  
  193. ;; Return the coefficient of the term of the polynomial e that
  194. ;; contains the highest power of x. When x = [x1,x2,...,xn], return 
  195. ;; lcoeff(lcoeff( ... (lcoeff(e,x1),x2),...,xn)...)
  196.  
  197. (defun $lcoeff (e &optional v)
  198.   (require-nonlist e "$lcoeff" 1)
  199.   (if (null v) (setq v (non-constant-ratvars e)))
  200.   (lcoeff (require-poly (myrat e) v "$lcoeff" 1) 
  201.       (require-list-of-symbols v "$lcoeff" 2)))
  202.   
  203. (defun lcoeff (e x)
  204.   (if (null x) e (lcoeff ($ratcoef e (car x) ($hipower e (car x))) (cdr x))))
  205.  
  206. ;; Return the coefficient of the term of the polynomial e that
  207. ;; contains the least power of x. When x = [x1,x2,...,xn], return 
  208. ;; lcoeff(lcoeff( ... (lcoeff(e,x1),x2),...,xn)...)
  209.  
  210. (defun $tcoeff (e &optional v)
  211.   (require-nonlist e "$tcoeff" 1)
  212.   (if (null v) (setq v (non-constant-ratvars e)))
  213.   (tcoeff (require-poly (myrat e) v "$tcoeff" 1) 
  214.       (require-list-of-symbols v "$tcoeff" 2)))
  215.   
  216. (defun tcoeff (e x)
  217.   (if (null x) e (tcoeff ($ratcoef e (car x) ($lowpower e (car x))) (cdr x))))
  218.  
  219. ;; Return the degree of the symbol x in the polynomial p.  When
  220. ;; x is a list, degree(p, [x1,x2,...,xn]) returns
  221. ;;       degree(p,x1) + degree(lcoeff(p, x1),[x2,...xn]).  
  222. ;; Finally, degree(p,[]) returns 0.
  223.  
  224. (defun $degree (p x)
  225.   (degree (require-poly (myrat p) x "$degree" 1) 
  226.       (require-list-of-symbols x "$degree" 2)))
  227.  
  228. (defun degree (p x)
  229.   (if (null x) 0
  230.     (add ($hipower p (car x)) (degree (lcoeff p `(,(car x))) (cdr x)))))
  231.  
  232. ;; Return the total degree of the polynomial.  Four cases:
  233. ;;   (a) total_degree(p) returns the total degree of the polynomial
  234. ;;       in the variables listofvars(p).
  235. ;;   (b) total_degree(p,x), where x isn't a list returns the 
  236. ;;       total_degree of p in the variable x.
  237. ;;   (c) total_degree(p,[x1,x2,...,xn]), where x = [x1,x2,...,xn] 
  238. ;;       returns the total_degree of p in the variables x1 thru xn.
  239. ;;   (d) total_degree(p,x1,x2,...xn), where the x's are symbols
  240. ;;       returns the total_degree of p in the variables x1 thru xn.
  241.  
  242. (defun $total_degree (p &optional v)
  243.   (if (null v) (setq v (non-constant-ratvars p)))
  244.   (setq v (require-list-of-symbols v "$total_degree" 2))
  245.   (total-degree (cadr (apply 'myrat `(,p ,@v)))))
  246.  
  247. (defun total-degree (e)
  248.   (cond ((consp (nth 2 e))
  249.      (+ (nth 1 e) (total-degree (nth 2 e))))
  250.     (t
  251.      (nth 1 e))))
  252.      
  253. (defun $lcm (p q)
  254.   (nth 1 ($divide (mul p q) ($gcd p q))))
  255.  
  256. ;; Compute the s-polynomial of f and g.  For a definition of the 
  257. ;; s-polynomial, see Davenport, Siret, and  Tournier, "Computer Algebra," 
  258. ;; 1988, page 100.
  259.  
  260. (defun $spoly (f g v)
  261.   (setq v (cons '(mlist simp) (require-list-of-symbols v "$spoly" 3)))
  262.   (setq f (myrat f))
  263.   (setq g (myrat g))
  264.   (let ((fp ($lterm f v))
  265.     (gp ($lterm g v)))
  266.     (mul ($lcm fp gp) (add (div f fp) (mul -1 (div g gp))))))
  267.  
  268. (defun $lterm (p &optional v)
  269.   (if (null v) (setq v (non-constant-ratvars p)))
  270.   (lterm (require-poly (myrat p) v "$lterm" 1) 
  271.      (require-list-of-symbols v "$lterm" 2)))
  272.  
  273. (defun lterm (p v)
  274.   (cond ((null v) p)
  275.     (t
  276.      (let* ((vo (car v))
  277.         (n ($hipower p vo)))
  278.        (lterm (mult ($ratcoef p vo n) (power vo n)) (cdr v))))))
  279.  
  280. (defun $get_exponents (p x)
  281.   (setq x (require-list-of-symbols x "$get_exponents" 2))
  282.   (let ((acc))
  283.     (setq p (myrat p))
  284.     (require-poly p (cons '(mlist simp) x) "$get_exponents" 1)
  285.     (dolist (xi x (cons '(mlist simp) (nreverse acc)))
  286.       (push ($hipower p xi) acc)
  287.       (setq p ($lcoeff p xi)))))
  288.  
  289. ;; Return true iff and only if e is a polynomial in the variables var.
  290.  
  291. (defun $polynomialp (e &optional vars cp ep)
  292.   (declare (ignore cp ep))
  293.   (if (null vars) (setq vars (non-constant-ratvars e)))
  294.   (setq vars (require-list-of-symbols vars "$polynomialp" 2))
  295.   (setq vars `((mlist simp) ,@vars))
  296.   (and (every #'(lambda (x) (or ($variablep x) ($ratfreeof vars x) 
  297.                 ($constantp x))) 
  298.           (cdr ($showratvars e)))
  299.        (not ($taylorp e)) ($ratfreeof vars ($ratdenom e))))
  300.  
  301.  
  302.      
  303.   
  304.   
  305.